De novo variants in the RNU4-2 snRNA cause a frequent neurodevelopmental syndrome

Around 60% of individuals with neurodevelopmental disorders (NDD) remain undiagnosed after comprehensive genetic testing, primarily of protein-coding genes1. Large genome-sequenced cohorts are improving our ability to discover new diagnoses in the non-coding genome. Here we identify the non-coding RNA RNU4-2 as a syndromic NDD gene. RNU4-2 encodes the U4 small nuclear RNA (snRNA), which is a critical component of the U4/U6.U5 tri-snRNP complex of the major spliceosome2. We identify an 18 base pair region of RNU4-2 mapping to two structural elements in the U4/U6 snRNA duplex (the T-loop and stem III) that is severely depleted of variation in the general population, but in which we identify heterozygous variants in 115 individuals with NDD. Most individuals (77.4%) have the same highly recurrent single base insertion (n.64_65insT). In 54 individuals in whom it could be determined, the de novo variants were all on the maternal allele. We demonstrate that RNU4-2 is highly expressed in the developing human brain, in contrast to RNU4-1 and other U4 homologues. Using RNA sequencing, we show how 5′ splice-site use is systematically disrupted in individuals with RNU4-2 variants, consistent with the known role of this region during spliceosome activation. Finally, we estimate that variants in this 18 base pair region explain 0.4% of individuals with NDD. This work underscores the importance of non-coding genes in rare disorders and will provide a diagnosis to thousands of individuals with NDD worldwide.


De novo variants in the RNU4-2 snRNA cause a frequent neurodevelopmental syndrome
Around 60% of individuals w it h n eu ro de ve lo pmental disorders (NDD) remain undiagnosed after comprehensive genetic testing, primarily of protein-coding genes 1 .Large genome-sequenced cohorts are improving our ability to discover new diagnoses in the non-coding genome.Here we identify the non-coding RNA RNU4-2 as a syndromic NDD gene.RNU4-2 encodes the U4 small nuclear RNA (snRNA), which is a critical component of the U4/U6.U5 tri-snRNP complex of the major spliceosome 2 .We identify an 18 base pair region of RNU4-2 mapping to two structural elements in the U4/U6 snRNA duplex (the T-loop and stem III) that is severely depleted of variation in the general population, but in which we identify heterozygous variants in 115 individuals with NDD.Most individuals (77.4%) have the same highly recurrent single base insertion (n.64_65insT).In 54 individuals in whom it could be determined, the de novo variants were all on the maternal allele.We demonstrate that RNU4-2 is highly expressed in the developing human brain, in contrast to RNU4-1 and other U4 homologues.Using RNA sequencing, we show how 5′ splice-site use is systematically disrupted in individuals with RNU4-2 variants, consistent with the known role of this region during spliceosome activation.Finally, we estimate that variants in this 18 base pair region explain 0.4% of individuals with NDD.This work underscores the importance of non-coding genes in rare disorders and will provide a diagnosis to thousands of individuals with NDD worldwide.
Despite increasingly powerful genomic and analytic approaches for the diagnosis of rare developmental disorders, around 60% of individuals with NDD remain without an identified genetic diagnosis after genomic testing with current methods 1 .So far, most known disease-causing variants are in the roughly 1.5% of the genome that directly encodes proteins 3 .By contrast, the non-coding genome (which makes up the remaining 98.5%) has been relatively unexplored, especially regions far from protein-coding genes.Large-scale systematic application of genome sequencing to clinical populations has increasingly enabled investigation of the contribution of variants in non-coding regions to genetic disorders 4 .
Non-coding RNAs, which comprise 37.4% of processed exonic RNA sequence in humans 5 , include important regulators of biological processes with diverse roles across cells and tissues 6 .snRNAs are a subcategory of non-coding RNAs that are key components of the spliceosome 7 .snRNAs complex with a multitude of proteins and other snRNA species in small nuclear ribonucleoprotein (snRNP) complexes to mediate the removal of introns from pre-messenger RNA transcripts 8 .Many spliceosome components have demonstrated roles in human disorders, including two snRNA components of the minor spliceosome: RNU12 variants cause autosomal recessive early-onset cerebellar ataxia 9 , whereas RNU4ATAC variants cause an autosomal recessive multisystem congenital disorder including microcephaly, growth retardation and developmental delay (eponyms include Taybi Linder 10 , Lowry-Wood 11 and Roifman syndromes 12 ).
Here we identify variants in RNU4-2, which encodes the U4 snRNA component of the major spliceosome, in an autosomal dominant disorder.Using a cohort of 8,841 probands with genetically undiagnosed NDD in the Genomics England 100,000 genomes project (GEL) 4 , we identify variants in a critical 18 base pair (bp) region in the centre of RNU4-2 associated with a severe neurodevelopmental phenotype and estimate that variants in this region explain around 0.4% of individuals with neurodevelopmental disorders (NDD).We demonstrate that variants in this region are severely depleted from large population datasets.We show that NDD variants map to critical structural elements in the U4/U6 complex that are important to correctly position U6 ACAGAGA to receive the 5′ splice site during initial spliceosome activation, and detail the expression of RNU4-2 through brain development.
The n.64_65insT variant is not found in 76,215 genome-sequenced individuals in gnomAD v.4.0 (ref.14), or in 245,400 individuals in the All of Us dataset 15 .It is seen in a single individual in the UK Biobank 16 (allele frequency of 1.02 × 10 −6 ) with a variant allele balance consistent with a true variant (23 reference and 18 (44%) alternate reads).This individual has an International Classification of Diseases 10th Revision (ICD-10) code for 'personal history of disease of the nervous system and sense organs' but no further phenotype data to assess a potential NDD diagnosis (Supplementary Table 1).
Given the high occurrence rate of this recurrent insertion, we wanted to rule out the possibility that it is a sequencing or mapping error, despite the overwhelming evidence of phenotype enrichment.Notably, the variant is a single A insertion after a run of four Ts, ruling out the most common cause of sequencing error for indels, polymerase slippage in homopolymer repeats.The variant calls were all high quality according to both analysis of quality metrics (Supplementary Fig. 1) and manual inspection on the Integrative Genomics Viewer (IGV) (Supplementary Fig. 2).The genomic region surrounding the insertion and RNU4-2 maps uniquely to a single region of the genome with short-read sequencing in GRCh38 and T2T CHM13v2.0/hs1.Finally, sequencing reads aligned to RNU4-2 map with good quality (average 96 reads with mapping quality scores greater than 20; Supplementary Fig. 3).

n.64_65insT is in a highly constrained region
The recurrent n.64_65insT variant resides within the central region of RNU4-2, towards the 5′ end of an 18 bp region that is depleted of variants in population datasets compared with the rest of the gene (26% of all possible single nucleotide variants (SNVs) observed in UK Biobank compared to a median of 78% across the rest of the gene; Fig. 1a and Extended Data Fig. 2a).On the basis of the population variant data, we defined a critical, highly constrained region as chr.12: 120291825-120291842.We refer to this as the 'critical region' throughout the rest of the paper.
We searched for variants across this region in GEL, and also in further cohorts containing undiagnosed individuals with NDD (Methods).In total, we identified 115 individuals with variants across this region, including 61 individuals in GEL (60 probands and one more sibling) and 54 from extra cohorts (Fig. 1b and Extended Data Table 1).For 86 of the 115 individuals, sequencing data for both parents were available to confirm that the variants had arisen de novo.Where possible, we used nearby variants to determine the parental allele of origin of the variants.For 54 individuals in whom this could be confidently resolved (46 with n.64_65insT; three with other insertion variants; five with SNVs), all 54 were present on the maternal allele.In one individual the n.65A>G variant appeared to be mosaic in the mother (54 reference and eight alternate reads) and in another an SNV was maternally inherited (n.76C>T).This analysis also enabled us to determine the likely de novo occurrence for five more individuals in whom only one parent was sequenced.Sanger sequencing was used to confirm the presence of the variant in eight individuals with the n.64_65insT variant.For seven of the eight, absence from both parents was also confirmed.In the eighth, the variant was confirmed as absent from the single available parent.In three families, the n.64_65insT variant was identified as occurring de novo in both short-and long-read trio sequencing.
Most of the 115 individuals have the initial n.64_65insT variant (n = 89; 77.4%).Five of the 11 extra variants are also single base insertions, including n.77_78insT (GRCh38:chr.12: 120291826:T:TA), which is seen in six individuals, two of whom are affected siblings.Single base insertion variants in this region are strongly enriched in individuals with NDD: 54 out of 8,841 (0.61%) GEL undiagnosed NDD probands (55 out of 10,388 individuals) have single base insertions compared to two out of 490,132 individuals in the UK Biobank (OR = 1,531; 95% confidence interval (CI): 404-16,384; Fisher's P = 3.3 × 10 −92 ).
Aside from insertions, there is also a modest enrichment of SNVs in GEL NDD probands across the critical region (undiagnosed NDD: six out of 8,841; UK Biobank: 35 out of 490,132; OR = 9.51; 95%CI: 3.27-22.8;Fisher's P = 8.16 × 10 −5 ).We identified 15 individuals across cohorts with SNVs in this region (Extended Data Table 1; 10 confirmed de novo), all with phenotypes consistent with individuals with insertion variants.The identified SNVs cluster with the two regions harbouring insertion variants at the extreme ends of the 18 bp critical region (Fig. 1).Conversely, SNVs in the central portion (particularly at nucleotides 71-74) are observed in both non-NDD individuals in GEL (n = 2) and population controls, although all at low frequencies (Extended Data Table 1).Across the remainder of RNU4-2 there is no significant enrichment of variants in undiagnosed NDD probands when compared with non-NDD probands (194 out of 7,519 undiagnosed NDD; 521 out of 19,428 non-NDD in GEL aggregated variant dataset 17 ; OR = 0.96; 95%CI: 0.81-1.14;Fisher's P = 0.67).
In total, we identify variants in this 18 bp region in 115 individuals with NDD.This includes 60 out of 8,841, or 0.68%, of all genetically undiagnosed NDD probands in GEL.By contrast, variants in this region are observed in 39 out of 490,132 (0.008%) individuals in the UK Biobank (OR = 85.8; 95%CI: 56.4-131.6;Fisher's P = 1.84 × 10 −78 ).As most individuals in GEL have had genetic testing before recruitment, we cannot estimate the overall prevalence of RNU4-2 variants in all cause NDD from this cohort.Instead, if we assume a diagnostic yield of 40% before defining our GEL undiagnosed NDD cohort, consistent with recent reports 1 , we estimate that variants in RNU4-2 could explain 0.4% of all NDD (calculated as 60 from an effective cohort size of 14,735 (8,841 × 1/0.6)).
U4 snRNA binds to U6 snRNA through extensive complementary base pairing in the U4/U6.U5 tri-snRNP complex of the major spliceosome.Unwinding of U6 from U4 is essential to generate the catalytically active spliceosome 2 .The 18 bp critical region in RNU4-2 maps to a region of U4 between the stem I region of complementary base pairing to U6 and the 3′ stem-loop structures (nucleotides 62-79; Fig. 1c).This region is known to be loaded into the active site of the SNRNP200-encoded BRR2 helicase, which mediates spliceosome activation by unwinding the U4/U6 duplex 2 .The highly recurrent n.64_65insT variant is within a previously described 'quasi pseudoknot', or T-loop, structure 18 (Fig. 1d).The region spanning nucleotides 76 to 78, where the recurrent n.77_78insT variant resides, is involved in base pairing with U6 in stem III (ref.19) (Fig. 1d).Both regions are thought to stabilize U4/U6 pairing and accurately position the U6 ACAGAGA sequence to receive the 5′ splice site during initial spliceosome assembly, whereas U4 nucleotides in stem III are part of the loading site for BRR2.Nearby regions that are predicted to have important roles, such as the U4/U6 stem I binding region, are not enriched for variants in NDD probands.

RNU4-2 variants disrupt 5′ splice-site use
Given the importance of U4 snRNA in the spliceosome and previous observations of global disruption to splicing observed in other spliceosomopathies 20 , we analysed RNA sequencing (RNA-seq) data from blood samples for five individuals from GEL.Three of these individuals have the highly recurrent n.64_65insT variant, another has the other recurrent insertion, n.77_78insT and the final individual has an SNV (n.78A>C).The five individuals with RNU4-2 variants had more abnormal splicing events than 378 control individuals with non-NDD phenotypes (mean 21.6 versus 4.5; Wilcoxon P = 0.0126), but this was not significant Article after correcting for multiple testing.There was no difference in the number of genes that were significant outliers for expression (mean 1.8 versus 5.7; Wilcoxon P = 0.94; Extended Data Table 2).
Consistent with the importance of the critical region in 5′ splice-site recognition, the most pronounced difference was observed for abnormal splicing events corresponding to increased use of unannotated 5′ splice sites (mean 8.8 events in individuals with RNU4-2 variants compared with 0.7 in both 378 unmatched controls and ten controls matched on genetic ancestry, sex and age at consent; Wilcoxon P = 4.0 × 10 −5 and P = 5.7 × 10 −3 , respectively; Fig. 2a and Extended Data Table 2).The individual with the SNV was not notably different from the four individuals with single base insertions (three significant events).Sequence motif analysis showed an increase in T at the +3 position and an increase in C at the +4 and +5 positions in the unannotated 5′ splice sites that were significantly increased in individuals with RNU4-2 variants compared to decreased canonical sites (Fig. 2c).These three positions of the 5′ splice site (+3, +4 and +5), which shift away from consensus in individuals with RNU4-2 variants, pair directly with the U6 ACAGAGA region during spliceosome activation (Fig. 2d).
Of all the detected abnormal splicing events, 12 of these were shared by two or more individuals with RNU4-2 variants (Extended Data Table 3).Eleven of these 12 events (91.6%) corresponded to an increase in unannotated 5′ splice-site use.None of these shared events were identified in any of the 378 controls.By contrast, when randomly sampling five control individuals across 10,000 permutations, the mean number of events shared by two or more individuals was 0.007, significantly fewer than the 12 in RNU4-2 individuals (permutation P < 1 × 10 −4 ; Fig. 2b).Five of the genes implicated in the 12 shared events are in the Developmental Disorders Genotype-to-Phenotype (DDG2P) database 21 and/or were associated with NDD in a previous    38 and are present in the RNA-seq reads from human prefrontal cortex in BrainVar.d, The structure of U4 and U6 snRNAs resolved by cryogenic electron microscopy 18 .U4 residues in the critical region are labelled with the reference nucleotide and numbered according to the position along the RNA (for example, U62 indicates a uracil residue in the reference sequence at position 62).Created using publicly accessible coordinates from the RCSB Protein Data Bank 39  large-scale analysis 22 (NDUFV1, H2AC6, JMJD1C, MAP4K4 and SF1).Ten of the 12 shared events affect the protein-coding sequence, with four predicted to cause a frameshift (Extended Data Table 3).Collectively, these results indicate a systematic shift in 5′ splice-site use in individuals with RNU4-2 variants compared to controls.Future work should assess these patterns in a more disease-relevant tissue (for example, brain) or in induced pluripotent stem cell-derived neuronal cells or organoid models.At present, RNA from more tissues from affected individuals is not available.

Characterizing the RNU4-2 NDD syndrome
To characterize the phenotypic spectrum associated with variants in RNU4-2, we collected detailed phenotypic information for a subset of 49 individuals (42 with n.64_65insT, three with other single base insertions, and four with SNVs; Table 1 and Supplementary Table 2).Using these data, we find the RNU4-2 syndromic NDD to be characterized by moderate to severe global developmental delay (four children with SNVs with moderate delay) and intellectual disability in all individuals.Most (83%) achieved ambulation but at a delayed age (average 3.4 years, range 17 months to 7.5 years) with some noted to have a wide-based or ataxic gait.Only three individuals (two with an SNV) had fluent speech, some had a few words and most were non-verbal.All but three were reported to have dysmorphic facial features.These facial features varied but consisted of a myopathic face with deep set eyes (some widely spaced and some narrowly spaced), epicanthus, wide nasal bridge, anteverted nares or underdeveloped ala nasi, large cupped ears (some posteriorly rotated), full cheeks, a distinctive mouth with full lips with downturned corners, high arched palate and a large or protruding tongue (Fig. 3).In comparison to the single base insertions, children with SNVs had fewer reports of severe global developmental delay (zero out of four versus 34 out of 40, Fisher's P = 0.0015).
Associated growth and neurodevelopmental phenotypes present in more than or equal to 75% of individuals include short stature, microcephaly (mostly congenital), speech abnormalities (mostly non-verbal), hypotonia and seizures.Seizures had variable onset with a median of 3 years and ranging from the first year of life to 10 years of age (spanning infantile spasms, focal seizures and generalized tonic-clonic seizures, febrile seizures and status epilepticus).Brain magnetic resonance imaging showed a spectrum of abnormalities in most individuals, most frequently reduced white matter volume, hypoplasia of the corpus callosum, ventriculomegaly, delayed myelination and other non-specific abnormalities of the white matter.Involvement of several organ systems was reported for all individuals (with fewer systems reported as involved in individuals with SNVs), often including visual (optic nerve hypoplasia, cortical blindness, strabismus, nystagmus), gastrointestinal (constipation, reflux, feeding issues with need for a gastrostomy tube, poor growth) and bone and/or skeletal abnormalities (osteopenia, recurrent fractures, scoliosis, kyphosis, hip dysplasia) and in fewer individuals, hearing, endocrine (hypothyroidism, growth hormone deficiency, pituitary abnormalities), limb, sleep, genitourinary, dental, cardiac and cutaneous concerns (Table 1 and Supplementary Table 2).No significant differences were noted in the presentation of male versus female individuals.

Exome sequencing rarely finds n.64_65insT
Most individuals with NDD who undergo genetic testing at present have exome rather than genome sequencing.Although RNU4-2 is not directly captured by exome sequencing panels, there is a chance that off-target reads may map to the 18 bp critical region of RNU4-2 and enable detection of variants in this region.To investigate this, we analysed individuals who are included in GEL and also have exome sequencing data in the Deciphering Developmental Disorders (DDD) cohort 1 .Across the DDD cohort, 3,408 out of 13,450 individuals (25.3%) have at least one read mapping to the position of the n.64_65insT variant (Extended Data Fig. 3).The maximum number of mapping reads in any individual was five, which is below standard thresholds used to identify heterozygous variants.Of 1,755 individuals in both GEL and DDD, 22 have the n.64_65insT variant (1.3%).Two of the 22 individuals (9.1%) each have a single read at the variant position in the exome sequencing data from DDD, but in each case it is identical to the reference sequence.The other 20 individuals have no reads mapping to RNU4-2.Nevertheless, others have reported successful identification and subsequent experimental validation of the n.64_65insT variant identified initially only on one or two sequencing reads (public communication on X/Twitter with S. Laurie and K. Platzer).These analyses suggest that although it is possible to identify individuals who may have variants in RNU4-2 through exome sequencing data, the sensitivity of this approach is very low.Any variants identified through this approach will also need independent confirmation.

Evaluating other snRNAs in NDD
Given the identified importance of RNU4-2 in NDD, we sought to determine whether other snRNA genes with no known association to NDD could also harbour new diagnoses.We investigated 28 snRNA genes that are expressed in the brain, using multiple approaches (Extended Data Table 4).First, we tested for overall enrichment of de novo variants in undiagnosed NDD probands compared to non-NDD probands across each snRNA with at least two identified de novo variants in probands with undiagnosed NDD (n = 14) using the high-confidence de novo callset in GEL.Of the 13 genes other than RNU4-2, none showed a significant enrichment of de novo variants in undiagnosed NDD probands (all Fisher's P > 0.15).Second, proposing that the burden of pathogenic variants in other snRNAs may be restricted to specific critical regions, as we see for RNU4-2, we used an 18 bp sliding window to identify snRNA regions that are depleted of variation in the UK Biobank compared to the overall variant burden across each gene.Notably, the regions with the highest depletion in RNU4ATAC correspond to two hotspots of pathogenic variants in ClinVar (chr.2: 121530923-121530946, chr.2: 121530984-121531007), however, the strength of the depletion in these regions is lower than in RNU4-2 (minimum normalized proportion of observed −0.11 and −0.2 versus −0.5 for the depleted region in RNU4-2), consistent with lower selection acting on variants in RNU4ATAC that cause recessive disorders.We identified 14 regions in 13 unique snRNAs with a deviation from the median number of SNVs across the full gene of at least 20% (Fig. 4 and Extended Data Table 5).We repeated our de novo variant enrichment test in regions with at least two de novo variants in undiagnosed NDD probands (n = 3).Only the conserved region in RNU4-2 was significant (Fisher's P = 9.31 × 10 −11 ; undiagnosed NDD probands n = 37, non-NDD probands n = 0; all other tests Fisher's P > 0.25).
Third, we looked for recurrent de novo variants in undiagnosed GEL NDD probands that were absent from diagnosed NDD probands, non-NDD probands and population controls.There are three de novo variants with an allele count greater than or equal to three in the GEL undiagnosed NDD cohort, two in RNU1-2 (chr.1: 16895992:C:T and chr.1: 16896002:A:G) and one in RNVU1-7 (chr.1: 148038767:G:A).However, all three variants are observed at comparable frequencies in non-NDD probands and are also found at relatively high frequencies in population controls (all variants' allele frequency greater than 0.5% in gnomAD v.4.0).
Finally, given that variants in RNU12 and RNU4ATAC are associated with recessive disease, we also tested for enrichment of homozygous and compound heterozygous variants in undiagnosed NDD probands compared to non-NDD probands.We observed a nominal enrichment of variants in RNU12 (11 probands with NDD versus two non-NDD probands; Fisher's P = 0.026), but this was not significant after correcting for multiple testing.We did not identify any significant associations across any other snRNA or when restricted to variants in our identified depleted regions (Extended Data Tables 4 and 5).

RNU4-2 is highly expressed in the brain
Humans have many genes that encode the U4 snRNA, although only two of these, RNU4-2 and RNU4-1, are highly expressed in the human brain (Supplementary Table 3).RNU4-2 and RNU4-1 are contiguous on chr.12, both 141 bp long and highly homologous, differing by four nucleotides (97.2% similarity).RNU4-1 has a similar depletion of variants in population cohorts in the centre of the RNA (Fig. 4), however, we do not observe an enrichment of variants in GEL in this central region (Extended Data Fig. 2b).There is a variant equivalent to our highly recurrent variant in RNU4-1 that is observed in six individuals in the UK Biobank dataset.There are no consistent phenotypes recorded in these six individuals (Supplementary Table 1).

Article
To investigate the reason for variants in RNU4-2, but not RNU4-1, causing NDD, we analysed the expression of both RNU4-1 and RNU4-2 in the brain.First, we analysed the expression patterns of both genes across many developmental stages using bulk RNA-seq data from 176 human prefrontal cortex samples in BrainVar 23 .The expression of RNU4-1 and RNU4-2 is tightly correlated (Supplementary Fig. 5), however, RNU4-2 is consistently expressed at a significantly higher level than RNU4-1 (Fig. 5a).Second, we assessed chromatin accessibility in the chromosome 12 locus containing both RNU4-1 and RNU4-2 using assay for transposase-accessible chromatin with sequencing (ATAC-seq) data from two human prenatal prefrontal cortex samples.These data show a notable chromatin accessibility signal around RNU4-2 and a much lower signal surrounding RNU4-1, again consistent with much higher expression of RNU4-2 in the brain (Fig. 5b).Overall, these data support the role of RNU4-2 as the main U4 transcript in the brain.

Discussion
Here we identified a highly constrained 18 bp region of RNU4-2 in which variants cause a severe neurodevelopmental phenotype.
We estimated that variants in this region could explain 0.4% of all NDD.As a comparison, the largest proportion of NDD explained by a single gene in 13,449 individuals in the DDD cohort 1 was 0.68% for ANKRD11, although we acknowledge that some genes and recognizable syndromes with longstanding associations (for example, MECP2, SCN1A, UBE3A) will be depleted from this cohort.The proportion of NDD explained by variants in RNU4-2 would be even higher if restricted to individuals with severe, syndromic NDD.This is consistent with the much lower rate of RNU4-2 variants in cohorts recruited primarily for autism spectrum disorder (for example, three out of 7,149; 0.042% across the Simons Simplex Collection (SSC) 24 , SPARK 25 and MSSNG 26 ).
Our findings underscore the value of large-scale genome sequencing datasets and the importance of considering variants outside protein-coding regions.This region, despite being within a highly conserved non-coding exon, is not included in commercially available clinical exome sequencing capture, which primarily targets protein-coding exons 5 .This discovery was empowered by the availability of increasingly large genome sequencing datasets from families affected by genetic disease around the world.Indeed, the scale and accessibility of the Genomics England dataset facilitated both the work reported here and a parallel discovery by an independent group 27 .The detailed phenotypic characterization included here will help prioritize individuals for targeted sequencing of RNU4-2.
For all individuals in whom we were able to confidently determine the parent of origin of the identified RNU4-2 variants (n = 54), the variants were observed to be on the maternal allele.This is in contrast to the well-established paternal bias observed for de novo small mutations 28 .The absence of any paternally derived variants in our cohort may be a consequence of negative selection in the male germline if RNU4-2 plays an important role during spermatogenesis.It may also be a consequence of imprinting, for example if variants on a highly expressed paternal allele are embryonic lethal, whereas those on a weakly expressed maternal allele are survivable but result in NDD.Further work is needed to test these hypotheses.
Most individuals in our cohort have the highly recurrent n.64_65insT variant.It is observed in 46 of 8,841 undiagnosed NDD probands in GEL.By contrast, the most recurrent protein-coding variant in a dataset of 31,058 individuals with developmental disorders 29 is observed in 36 individuals (0.12%; GRCh38:chr.11: 66211206:C:T; PACS1:p.Arg203Trp).The reasons for this high recurrence are unclear, but it could be driven by either a high endogenous mutation rate or positive selection in the germline.The latter has previously been described for so-called selfish mutations associated with paternal age effects 30 .One hypothesis is that germline selection is acting to increase the apparent frequency of the n.64_65insT variant, for example through meiotic drive effects or by accelerating oocyte maturation 31 .We do not see an association with maternal age for individuals with n.64_65insT in GEL (mean 30.2 compared to 29.7 across other NDD probands; Extended Data Fig. 4).
Alternatively, recurrence may be driven by a high mutation rate.This is consistent with the observed open chromatin state and very high expression of RNU4-2 (Fig. 5), as high levels of transcription are known to be correlated with increased mutation rate 32 .Hypermutability of short non-coding RNA genes, including snRNAs, has previously been documented 33,34 .Consistent with this, a high variant density is observed across RNU4-2 in the UK Biobank (Extended Data Fig. 5).Despite the high number of variants in RNU4-2 in the UK Biobank, there are no individuals with homozygous variants and all observed variants are very rare (maximum allele frequency 0.025%), consistent with strong negative selection acting on variants across RNU4-2.A high overall mutational burden does not, however, explain the high recurrence of this specific single base insertion.Local formation of secondary structure and base stacking is a known driver of biased small insertion mutations 35 .The high propensity of this region to form secondary structure when single-stranded may drive creation of this a, Levels of RNU4-1 (grey) and RNU4-2 (teal) expression in human dorsolateral prefrontal cortex at different developmental stages from BrainVar 23 .Coloured lines correspond to the Loess smoothed average with the shaded regions representing 95% CIs.Developmental stages are labelled with periods (1 to 12), spanning from embryonic development to late adulthood, that were defined previously 42 .b, ATAC-seq data from human prenatal prefrontal cortex shows substantially higher peaks of chromatin accessibility around RNU4-2 than RNU4-1.Data for both 18 and 19 gestational weeks (GW) are shown to demonstrate replication.
specific insertion.It is also possible that this variant is more compatible with live birth relative to other comparably recurrent mutations in the critical region.The n.64_65insT variant is one of six single base insertions that we observe in the 18 bp critical region in individuals with NDD, in a total of 100 individuals across cohorts.By contrast, single base insertions are very rare in population cohorts.Although we do also observe some SNVs in this region in individuals with NDD, our initial data suggest these SNVs may result in a milder phenotype.However, given this observation is based on only four fully phenotyped individuals, it needs to be confirmed in larger cohorts.Saturation mutagenesis experiments that test the impact of different length insertions and deletions as well as SNVs across the length of RNU4-2 will be important to understand the spectrum of deleterious mutations.The high proportion of single base insertion variants in individuals with NDD may indicate that steric conformational changes are needed to disrupt RNU4-2 function.Specifically, insertion of a single base into the T-loop or stem III regions may destabilize the U4/U6 interaction and/or alter the positioning of the U6 ACAGAGA sequence and potentially disrupt the correct loading of the 5′ splice site into the fully assembled spliceosome.This proposed effect is supported by the observed systematic disruption to 5′ splice-site use observed in RNA-seq data from five individuals with RNU4-2 variants.In particular, our observation that the +3, +4 and +5 positions of the 5′ splice site, which directly pair with the U6 ACAGAGA sequence, shift away from consensus at sites with increased use in individuals with RNU4-2 variants provides functional evidence that these variants disrupt accurate splice-site recognition during spliceosome activation.Further, variants in U6 snRNA and protein components of the spliceosome situated in the proximity of our RNU4-2 variants have recently been shown to alter 5′-splice-site selection by changing the preference for nucleotides that pair with the U6 snRNA ACAGAGA, consistent with this region being involved in subtle regulation of alternative splicing 36,37 .
Whereas two other snRNA genes, RNU12 and RNU4ATAC, have been linked to different phenotypes, both are components of the minor spliceosome and are associated with recessive disorders.By contrast, here we implicate variants in a major spliceosome snRNA in a dominant disorder.We further explored whether other snRNA genes could explain undiagnosed cases.We did not find any other snRNAs, or constrained subregions of snRNAs, that were significantly enriched for either de novo variants or recessively inherited variants in NDD cases when compared with non-NDD probands.We note, however, that these tests have low power given the small size of the genes and regions (mean 139.5 and 28.1 bp, respectively).Variants in the regions we delineated should also be investigated in other disease cohorts.
In summary, we identify RNU4-2 as a syndromic NDD gene, explaining roughly 0.4% of all individuals with NDD.Including RNU4-2 in standard clinical workflows will end the diagnostic odyssey for thousands of patients with NDD worldwide, and knowledge of the gene responsible for disease will enable investigation of potential treatments for these individuals.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-024-07773-7.
For most of our analyses, we used two previously defined datasets within GEL.First, a high-confidence set of de novo variants from 13,949 trios 13 .As of 13 March 2024, this dataset includes 12,554 probands with consent: 5,426 probands with undiagnosed NDD, 2,352 with diagnosed NDD and 4,776 non-NDD probands.De novo variants were filtered to those that pass the stringent_filter.Second, an aggregated variant call set (aggV2) 17 that contains 29,850 probands: 7,519 undiagnosed NDD, 2,903 diagnosed NDD and 19,428 non-NDD.

Expanded NDD cohort and clinical data collection
Clinical data were collected from research participants after obtaining written informed consent from their parents or legal guardians, with the study approved by the local regulatory authority.Samples were collected largely through personal communications (N.W., A.O.'D.-L., D.G.M.), as variants in this gene have not been prioritized in analysis.On entry into Matchmaker Exchange using the seqr 45  Clinical data were collected and summarized for features seen across the cohort.Written consent was obtained to publish all photographs included in Fig. 3.

Generating 1,000 random intergenic sequences
Using the bedtools (v.2.31.0) subtractBed function 47 we retrieved regions on chromosome 12 that do not overlap with RefSeq transcripts aligned by the National Center for Biotechnology Information.
We further removed regions within 10 kbp of an annotated transcript and restricted the remaining regions to those at least 141 bp in length (n = 611).We further removed regions overlapping the centromere.We then generated a set of 1,000 random sequences from each intergenic region and then randomly selected 1,000 non-overlapping regions from these.

Identifying human snRNA genes
We extracted genes with snRNA biotypes from Ensembl genome annotation v.111.We filtered out known pseudogenes (that is, with gene names marked with 'P' or identified through manual curation).For each remaining gene, we used BrainVar 23 RNA-seq expression data to calculate the mean counts per million (CPM) value across the gene.We selected only genes with mean CPM value across all BrainVar samples greater than five, resulting in a dataset of 28 snRNA genes.

Assessing variant depletion
Given the high mutability of RNU4-2 and other snRNA genes, coupled with strong selection pressures on variants, we did not think that conventional mutational models would be well calibrated to assess variant depletion.Instead, we devised a sliding window-based strategy to identify regions within snRNA genes that are relatively depleted of SNVs.We split genes into 18 bp sliding windows (chosen as it is the size of the region defined in RNU4-2) and tallied the number of SNVs observed in UK Biobank 500k genome sequencing data within that window, divided by the total number of possible SNVs (that is, 18 × 3).The proportion of possible SNVs observed in each window was normalized to the median across all sliding windows in that gene (that is, the per-gene median proportion observed was subtracted from each value).Depleted regions were defined as those spanning windows with a deviation from the per-gene median of at least 20%, that is, normalized observed proportion of possible SNVs less than −0.2.The same calculation was performed on 1,000 randomly selected 141 bp intergenic regions on chr. 12 (above).A one-way approximative (Monte Carlo) Fisher-Pitman test was conducted to show the median observed proportion of possible SNVs was significantly higher for RNU4-1 and RNU4-2 compared to the distribution in the 1,000 random regions.

RNA-seq of individuals with RNU4-2 variants
Blood was collected from a subset of 100,000 Genomes Project probands in PaxGene tubes to preserve RNA at the time of recruitment.RNA was extracted, depleted of globin and ribosomal RNAs, and subjected to sequencing by Illumina using 100 bp paired-end reads, with a mean of 102 million mapped reads per individual.Alignment was performed using Illumina's DRAGEN pipeline (v.3.8.4).FRASER2 (ref.40) and OUTRIDER 48 were used to detect abnormal splicing events and expression differences with samples run in batches of 500, both run using the DROP pipeline 49 (v.1.3.3).Significant outlier events were identified as those with a false discovery rate adjusted P < 0.05.The number of outlier events detected in five individuals with RNU4-2 variants was compared to two different control sets: (1) ten individuals matched (two per RNU4-2 individual) on genetic ancestry, sex and approximate age at consent and that did not have any NDD phenotypes; (2) 378 individuals with more than 60 million mapped reads, age below 18 years and with no NDD phenotypes.Sequence logo plots in Fig. 2 and Supplementary Fig. 6 were created in R (v.4.0.2) using the ggseqlogo package.

Assessing the sensitivity to detect the n.64_65insT variant in exome sequencing data
We used a Python script that uses samtools mpileup to retrieve the coverage and base change at the n.64_65 critical locus to identify putative carriers of the insertion (https://github.com/francois-lecoquierre/genomics_shortcuts/blob/main/find_RNU4-2_recurrent_variant.py).This was applied to exome sequencing data (32,681 CRAM files from probands and parents) from the DDD cohort.

Analysing RNU4-2 and RNU4-1 expression
We used the BrainVar 23 dataset to assess patterns of whole-gene expression of RNU4-2 and RNU4-1 in the human cortex across prenatal and postnatal development.This dataset includes bulk-tissue RNA-seq data from 176 de-identified postmortem samples of the dorsolateral prefrontal cortex (DLPFC, n = 167 older than ten postconception weeks) or frontal cerebral wall (n = 9 younger than ten postconception weeks), ranging from six postconception weeks to 20 years of age.The 100 bp paired-read RNA-seq data from BrainVar were aligned to the GRCh38.p12human genome using STAR aligner 50 (v.2.4.2a), and gene-level read counts for GENCODE v.31 human gene definitions were calculated with DEXSeq 51 (v.1.50.0) and normalized to CPM 52 .

Prenatal prefrontal cortex ATAC-seq data
Methods of generating ATAC-seq have been described previously 53 , which is the source of the data shown here.Briefly, fresh prenatal (18 and 19 gestational weeks) brain samples were dissected within 2 h of elective termination to extract the entire telencephalic wall, from the ventricular zone to the meninges.Intact nuclei were isolated by manually douncing the tissue on ice using a loose pestle douncer then lysed on ice for 10 min by adding a solution with 0.1% NP-40.Nuclei were spun down by centrifugation then resuspended and exposed to Tagmentation Enzyme for 30 min at 37 °C.The ATAC-seq library was generated using Illumina barcode oligos, amplified by high-fidelity PCR and sequenced on the Illumina HiSeq 2500 using paired-end sequencing.Reads were aligned to GRCh38 using the ENCODE ATAC-seq pipeline with default parameters 54 .A UCSC Browser track of per nucleotide ATAC-seq counts was used to assess the region around RNU4-2 and RNU4-1.

Burden testing and statistical analysis
The enrichment of both de novo variants and homozygous and/or compound heterozygous variants across each of 28 snRNA genes and 14 constrained subregions was assessed in undiagnosed NDD probands compared to non-NDD probands.De novo variants were identified from GEL's high-confidence de novo callset.Homozygous and compound heterozygous variants were identified from genotyping data in individual participants' VCF files.Homozygous variants were identified as variants that are heterozygous in both parents and homozygous in offspring.To identify compound heterozygous variants in a gene or region, we assessed whether there are greater than or equal to 1 paternally inherited heterozygous variant and greater than or equal to one maternally inherited heterozygous variant in the offspring.Multiallelic sites were excluded from this analysis.Homozygous variants and compound heterozygous variants were grouped together for burden testing.ORs and associated P values were calculated using a two-sided Fisher's exact test in R. A P value threshold of 0.0014 was used to assess statistical significance as a Bonferroni correction accounting for 35 tests.

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.A median of 76% of all possible SNVs in RNU4-2 are observed compared with 13% on average in the intergenic sequences of the same length (141 bp; P < 0.001, Monte Carlo Fisher-Pitman test).

Extended Data Table 4 | Burden testing across snRNAs
Genomic coordinates of, and burden testing results for snRNA genes in 5,426 undiagnosed NDD probands against 4,776 non-NDD probands.

April 2023
Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies.Here, indicate whether each material, system or method listed is relevant to your study.If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Fig. 1 |
Fig. 1 | A highly structured 18 bp region of RNU4-2 that is critical for BRR2 helicase activity is enriched for variants in NDD and depleted in population cohorts.a, Allele counts of de novo variants in 8,841 undiagnosed NDD probands in GEL (top; teal) and the UK Biobank cohort (bottom; grey) across RNU4-2.The 18 bp critical region, which is depleted of variants in the UK Biobank, is marked by a horizontal bar at the top of the plot.b, Allele counts of further variants identified in individuals with NDD in the critical 18 bp region.This includes 16 individuals with seven variants without sequencing data for both parents in GEL and variants identified in individuals from the following extra cohorts (Methods): NHS GMS (n = 19); MSSNG (n = 2); SSC (n = 1); GREGoR (n = 10) and UDN (n = 6); from personal communication or Matchmaker Exchange (n = 16).c, Schematic of U4 (teal) binding to U6 snRNA (grey).The 18 bp critical region is underlined.Nucleotides 142 to 145 of U4 (in blue) are not within the GENCODE transcript of RNU4-2 but are included in previous figures of the U4/U6 duplex in the literature on which this depiction is based38 and are present in the RNA-seq reads from human prefrontal cortex in BrainVar.d, The structure of U4 and U6 snRNAs resolved by cryogenic electron microscopy18 .U4 residues in the critical region are labelled with the reference nucleotide and numbered according to the position along the RNA (for example, U62 indicates a uracil residue in the reference sequence at position 62).Created using publicly accessible coordinates from the RCSB Protein Data Bank 39 (PDB structure 6QW6).In both c and d, single base insertions identified in individuals with NDD are shown by black arrows and positions of SNVs by orange nucleotides.
(PDB structure 6QW6).In both c and d, single base insertions identified in individuals with NDD are shown by black arrows and positions of SNVs by orange nucleotides.

Fig. 2 |
Fig. 2 | Individuals with RNU4-2 variants have systematic changes in 5′ splicesite use.a, Boxplots of the number of abnormal splicing events (detected by FRASER2, ref. 40) at unannotated 5′ splice sites.The individuals with RNU4-2 variants (n = 5 individuals) have significantly more outlier events than both controls with non-NDD phenotypes (n = 378 individuals) and controls matched on genetic ancestry, sex and age at consent (n = 10 individuals, two per case; Wilcoxon P = 4.0 × 10 −5 (W test statistic, 1,766) and P = 5.7 × 10 −3 (W test statistic, 45.5), respectively).Centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; maxima and minima represented as points.FDR, false discovery rate.b, The distribution of the number of abnormal splicing events at unannotated 5′ splice sites shared between two or more out of five randomly selected control individuals over 10,000 permutations (grey histogram).The number of shared events in individuals with RNU4-2 variants is indicated as a dotted teal vertical line (n = 11).c, DNA sequence motifs around 5′ splice sites with increased and decreased use in individuals with RNU4-2 variants.Each plot shows the proportion of sites with each base at each position.5′ splice sites with increased use (top) have an increase in T at the +3 position (eight out of 19 versus zero out of 36; Fisher's P = 6.2 × 10 −5 ; OR = Inf; 95%CI: 5.92-Inf) and an increase in C at the +4 (four out of 19 versus zero out of 36; Fisher's P = 0.011; OR = Inf; 95%CI: 1.37-Inf) and +5 (six out of 19 versus 1/36; Fisher's P = 0.0051; OR = 15.3;95%CI: 2.09-Inf)) positions compared to decreased 5′ splice sites (bottom).The consensus sequence at 5′ splice sites in matched annotation from NCBI and EMBL-EBI (MANE) transcripts 41 is shown in Supplementary Fig. 4. d, The structure of the U6 snRNA paired with the 5′ splice site after 5′ splice-site transfer.The three bases of the U6 ACAGAGA that directly pair with the 5′ splice site are shown in pink.The paired positions of the 5′ splice site (5′SS) are shown in green (A + 3 and A + 4) and yellow (G + 5).Statistical tests in a and c are one-sided with unadjusted P values.

Fig. 3 |Fig. 4 |
Fig. 3 | Clinical photographs showing facial features of affected individuals with variants in RNU4-2.All individuals shown have the n.64_65insT variant, except for individual 44 in o (n.68_69insA), individual 45 in p (n.64_65insG) and individual 48 in q (n.76C>T).a, Individual 1 at 12 years old.b, Individual 4 at 9 years old.c, Individual 7 at 13 years old.d, Individual 15 at 8 years old.e, Individual 21 at 3.5 years old.f, Individual 22 at 8 years old.g, Individual 23 at 13 years old.h, Individual 28 at 5 years old (left) and 9 years old (right).i, Individual 32 at 3 years old (left) and 12 years old (right).j, Individual 36 at 11 months old (left) and 8 years old (right).k, Individual 37 at 22 months old (left) and 16 years old (right).l, Individual 38 at 2.5 years old (left) and 10 years old (right).m, Individual 39 at 2 years old (left) and 12 years old (right).n, Individual 43 at 8 years old (left) and 12 years old (right).o, Individual 44 at 6 years old (left) and 19 years old (right).p, Individual 45 at 6 years old (left and centre) and 27 years old (right).q, Individual 48 at 22 months old.

Fig. 5 |
Fig.5| RNU4-2 is more highly expressed than RNU4-1 in the prefrontal cortex.a, Levels of RNU4-1 (grey) and RNU4-2 (teal) expression in human dorsolateral prefrontal cortex at different developmental stages from BrainVar23 .Coloured lines correspond to the Loess smoothed average with the shaded regions representing 95% CIs.Developmental stages are labelled with periods (1 to 12), spanning from embryonic development to late adulthood, that were defined previously42 .b, ATAC-seq data from human prenatal prefrontal cortex shows substantially higher peaks of chromatin accessibility around RNU4-2 than RNU4-1.Data for both 18 and 19 gestational weeks (GW) are shown to demonstrate replication.
node, one match was made with GeneMatcher 46 (C.D.).N.W. reviewed the National Health Service Genome Medicine Service (NHS GMS; v.3) dataset.Samples from NHS GMS were manually checked to remove duplicates with GEL.A.O.'D.-L.and S.S. reviewed the Broad Centre for Mendelian Genomics and the GREGoR consortium datasets and contacted the Undiagnosed Disease Network (UDN) through M.T.W. D.G.M. contacted extra local collaborators.Clinical collaboration requests were submitted to GEL to contact recruiting clinicans and collect extra phenotypic information.

Extended Data Fig. 1 |
HPO terms for individuals in GEL.(a) The proportion of individuals with human phenotype ontology (HPO) terms corresponding to phenotypes observed in ≥ 5 individuals with the n.64_65insT variant compared to all other individuals with NDD.Multiple HPO terms are significantly enriched in individuals with the n.64_65insT variant after Bonferroni adjustment (marked with a *) indicating that individuals with the n.64_65insT variant have more phenotypic similarity than the GEL NDD cohort as a whole.Multiple terms relating to global developmental delay, intellectual disability, hypotonia, seizure, microcephaly, autism, and short stature have been collapsed into single phenotypes.Of note, this figure relates only to HPO terms entered for each individual into GEL, which may be incomplete.Error bars indicate ±1 standard error.(b) Data plotted in panel (a) including statistics from twosided Fisher's exact tests.A P-value threshold of 2.94 × 10 −3 was used to assess statistical significance (Bonferroni adjusted for 17 tests).Extended Data Fig. 2 | Depletion of variants in the population in RNU4-2 and RNU4-1.(a; top) Distance to the median proportion of all possible SNVs that are observed in the UK Biobank in 18 bp sliding windows across the length of RNU4-2.A clear region of depletion compared to the rest of the gene is observed in the centre.(bottom) Log transformation of the mean Roulette 33 mutability across the 3 possible SNVs within a site.(b) Total allele frequency at each site of RNU4-1 in undiagnosed NDD probands in GEL (teal) and the UK Biobank cohort (grey).In contrast to RNU4-2, variants in RNU4-1 have higher allele frequencies.A similar region of depletion is seen in the centre of RNU4-1 (quantified in Fig. 4), but this is not enriched for variants in GEL NDD or non-NDD individuals.Extended Data Fig. 3 | Sequencing coverage in exome sequencing data.The number of sequencing reads covering the position of the n.64_65insT variant in 13,450 probands with exome sequencing in the DDD cohort.3,408/13,450 probands (25.3%) have at least one read at the position.Extended Data Fig. 4 | Comparison of parental age.Comparison of (a) paternal age for probands with fathers and (b) maternal age for probands with mothers recruited into GEL for individuals with variants in RNU4-2 (teal) and all other NDD probands (grey).Centre line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range.Individual data points, including outliers, are not shown due to Genomics England restrictions.NS: not significant.Paternal age: mean 33.1 vs 33.4 in individuals with RNU4-2 variants and other NDD probands, respectively (two-sided t-test P-value = 0.771; t = −0.29 (−2.41 -1.80)).Maternal age: mean 30.2 vs 29.7 in individuals with RNU4-2 variants and other NDD probands, respectively (two-sided t-test P-value = 0.505; t = −0.67 (−1.07 -2.15)).Extended Data Fig. 5 | Assessing variant density in the UK Biobank.Median proportion of possible SNVs observed in UK Biobank per 18 bp window across 1,000 intergenic regions on chromosome 12 (grey) and RNU4-1, RNU4-2 (teal).

Table 1 | Clinical features of 49 individuals with RNU4-2 variants Clinical feature
F, female; M, male; IUGR, intrauterine growth restriction; GDD, global developmental delay; ID, intellectual disability; ASD, autism spectrum disorder; MRI, magnetic resonance imaging; GORD, gastro-oesophageal reflux disease; G-tube, gastrostomy tube.a Maternal age available for 43 out of 49 individuals.b Paternal age available for 41 out of 49 individuals.c Denominator indicates the number of individuals for whom data were available.d One individual has a dual diagnosis in GJB2 that would account for the hearing loss.